Overview

Study Design

Research Question: Do damselfish (Dascyllus flavicaudus) enhance coral skeletal growth rates in Pocillopora spp., and does this effect vary with wound size?

Experimental Design:

Factors: Fish presence (Present/Absent) × Wound size (None/Small/Large)
Sample size: 72 coral fragments
Duration: 21 days
Location: Mo’orea, French Polynesia
Response variables:
- Buoyant weight (skeletal mass)
- Allometrically-scaled growth
- Growth rate per surface area per day

Analysis Overview

This document contains:

Data Import & Preparation
Allometric Scaling Analysis (size-corrected growth)
Surface Area Calibration
Growth Rate Models (mixed-effects)
Publication Figures
Comprehensive Diagnostics

Key Methodology: We use allometric scaling to correct for initial coral size differences, allowing fair comparison of growth rates across all fragments.

1. Setup & Data Import

Load Required Packages

# ==================== Core Packages ====================
library(tidyverse)      # Data wrangling and visualization
library(here)           # Relative file paths
library(janitor)        # Data cleaning

# ==================== Statistical Modeling ====================
library(lme4)           # Linear mixed-effects models
library(DHARMa)         # Residual diagnostics
library(broom.mixed)    # Tidy model outputs

# ==================== Visualization ====================
library(ggpubr)         # Publication-ready plots
library(gt)             # Publication-quality tables

cat("✓ Packages loaded successfully\n")

## ✓ Packages loaded successfully

Define Aesthetics

# Consistent color schemes for all figures
wound_palette <- c(
  "Large" = "#FF8B00",     # Orange
  "Small" = "#be8333",     # Tan  
  "No Wound" = "#485900"   # Olive green
)

fish_palette <- c(
  "Fish" = "#2b4b52",      # Teal
  "No Fish" = "#ad301a"    # Burnt red
)

Import Data

# Import buoyant weight data (skeletal mass measurements)
bw_raw <- read_csv(here("data", "fish_regen_buoyantweight (1).csv")) %>%
  clean_names() %>%
  # Remove outliers identified in QC
  filter(coral_id != 45) %>%
  filter(coral_id != 91)

# Import surface area data
sa <- read_csv(here("data", "fish_regen_mastersheet_wound closure_necrosis_sa - mastersheet.csv")) %>%
  clean_names() %>%
  mutate(coral_id = as.character(coral_id))

# Import tank metadata
tank_df <- read_csv(here("data", "Copy of fish_regen_mastersheet_wound closure_necrosis - mastersheet.csv")) %>%
  clean_names() %>%
  select(coral_id, tank) %>%
  mutate(coral_id = as.character(coral_id))

cat("✓ Data imported successfully\n")

## ✓ Data imported successfully

cat("  - Buoyant weight data:", nrow(bw_raw), "observations\n")

##   - Buoyant weight data: 134 observations

cat("  - Surface area data:", nrow(sa), "corals\n")

##   - Surface area data: 72 corals

cat("  - Tank assignments:", nrow(tank_df), "corals\n")

##   - Tank assignments: 72 corals

2. Data Preparation

Reshape Buoyant Weight Data

# Split data into initial and final measurements
bw_initial <- bw_raw %>%
  filter(date == "initial") %>%
  select(coral_id, wound, fish, initial_weight = bouyantweight_g)

bw_final <- bw_raw %>%
  filter(date == "final") %>%
  select(coral_id, final_weight = bouyantweight_g)

# Join into wide format (one row per coral)
bw_wide <- left_join(bw_initial, bw_final, by = "coral_id") %>%
  mutate(
    delta_mass = final_weight - initial_weight,
    coral_id = as.character(coral_id)
  )

cat("✓ Data reshaped to wide format\n")

## ✓ Data reshaped to wide format

cat("  - Initial weights:", sum(!is.na(bw_wide$initial_weight)), "corals\n")

##   - Initial weights: 67 corals

cat("  - Final weights:", sum(!is.na(bw_wide$final_weight)), "corals\n")

##   - Final weights: 67 corals

Merge with Surface Area Data

# Adjust surface area to account for wound area removed
sa_adjusted <- sa %>%
  select(-fish, -date_collected, -date_fragment) %>%
  mutate(
    # Subtract wound area: 3 cm² for large wounds, 1 cm² for small
    sa_col = ifelse(wound == "Large", sa_cal - 3, 
                    ifelse(wound == "Small", sa_cal - 1, sa_cal))
  ) %>%
  select(-wound)  # Remove to avoid duplication in join

# Merge buoyant weight with adjusted surface area
bw_sa_wide <- left_join(bw_wide, sa_adjusted, by = "coral_id") %>%
  mutate(
    wound = factor(wound, levels = c("No Wound", "Small", "Large")),
    tank = as.factor(tank)
  ) %>%
  mutate(tank = factor(tank, levels = sort(unique(as.numeric(as.character(tank))))))

cat("✓ Surface area data merged\n")

## ✓ Surface area data merged

cat("  - Mean initial SA:", round(mean(bw_sa_wide$sa_cal, na.rm = TRUE), 2), "cm²\n")

##   - Mean initial SA: 69.38 cm²

Merge with Tank Metadata

bw_merged <- left_join(bw_wide, tank_df, by = "coral_id") %>%
  mutate(
    wound = factor(wound, levels = c("No Wound", "Small", "Large")),
    tank = as.factor(tank)
  )

cat("✓ Tank metadata merged\n")

## ✓ Tank metadata merged

cat("  - Number of tanks:", nlevels(bw_merged$tank), "\n")

##   - Number of tanks: 6

3. Allometric Scaling Analysis

Theory: Why Allometric Scaling?

Problem: Larger corals grow faster simply because they’re bigger, not necessarily because they’re healthier.

Solution: Allometric scaling adjusts growth for initial size using the relationship:

\[\log(\text{Final Mass}) = \alpha + b \times \log(\text{Initial Mass})\]

Where \(b\) is the allometric coefficient. We then calculate size-corrected growth as:

\[\text{Scaled Growth} = \log\left(\frac{\text{Final Mass}}{\text{Initial Mass}^b}\right)\]

Prepare Data for Allometric Model

# Create dataset with log-transformed weights
log_df <- bw_merged %>%
  mutate(
    log_i_weight = log(initial_weight),
    log_f_weight = log(final_weight),
    growth_ratio = final_weight / initial_weight
  ) %>%
  drop_na(log_i_weight, log_f_weight, tank)

cat("✓ Allometric dataset prepared\n")

## ✓ Allometric dataset prepared

cat("  - Valid observations:", nrow(log_df), "corals\n")

##   - Valid observations: 67 corals

Fit Allometric Model

# Linear mixed-effects model: final ~ initial + random(tank)
model_lmer <- lmer(log_f_weight ~ log_i_weight + (1 | tank), data = log_df)

cat("=== Allometric Model Summary ===\n")

## === Allometric Model Summary ===

summary(model_lmer)

## Linear mixed model fit by REML ['lmerMod']
## Formula: log_f_weight ~ log_i_weight + (1 | tank)
##    Data: log_df
## 
## REML criterion at convergence: -515.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1853 -0.5898  0.0751  0.6736  2.0313 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  tank     (Intercept) 3.518e-07 0.0005931
##  Residual             1.899e-05 0.0043574
## Number of obs: 67, groups:  tank, 6
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.044640   0.005354   8.338
## log_i_weight 0.990376   0.001961 504.931
## 
## Correlation of Fixed Effects:
##             (Intr)
## log_i_weght -0.994

# Extract allometric coefficient (b-value) with 95% CI
b_value <- tidy(model_lmer, effects = "fixed", conf.int = TRUE) %>%
  filter(term == "log_i_weight")

cat("\n=== Allometric Coefficient (b-value) ===\n")

## 
## === Allometric Coefficient (b-value) ===

print(b_value)

## # A tibble: 1 × 7
##   effect term         estimate std.error statistic conf.low conf.high
##   <chr>  <chr>           <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 fixed  log_i_weight    0.990   0.00196      505.    0.987     0.994

Visualize Allometric Relationship

# Generate predicted values
fit_df <- log_df %>%
  select(log_i_weight) %>%
  distinct() %>%
  arrange(log_i_weight) %>%
  mutate(predicted = predict(model_lmer, newdata = ., re.form = NA))

# Create allometric scaling plot
allometric_plot <- ggplot(log_df, aes(x = log_i_weight, y = log_f_weight)) +
  geom_point(size = 2.2, shape = 21, fill = "steelblue", color = "black", 
             alpha = 0.8, stroke = 0.3) +
  geom_line(data = fit_df, aes(x = log_i_weight, y = predicted), 
            color = "darkred", linewidth = 1.2) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", 
              color = "gray60", linewidth = 0.8) +
  labs(
    title = "Allometric Scaling of Coral Buoyant Weight",
    subtitle = "Mixed-effects model with tank as random intercept",
    x = "Log Initial Weight (g)",
    y = "Log Final Weight (g)"
  ) +
  theme_pubr(base_size = 14)

print(allometric_plot)

# Save figure
ggsave(
  filename = here("figures", "growth", "allometric_scaling_plot.pdf"),
  plot = allometric_plot,
  width = 8, height = 6, dpi = 600, device = cairo_pdf
)

ggsave(
  filename = here("figures", "growth", "allometric_scaling_plot.png"),
  plot = allometric_plot,
  width = 10, height = 7.5, dpi = 300, bg = "white"
)

cat("✓ Allometric plot saved\n")

## ✓ Allometric plot saved

Display b-value Table

# Create formatted table of allometric coefficient
b_value_gt <- tidy(model_lmer, effects = "fixed", conf.int = TRUE) %>%
  filter(term == "log_i_weight") %>%
  transmute(
    Term = "Allometric slope (b-value)",
    Estimate = round(estimate, 3),
    `95% CI` = paste0("[", round(conf.low, 3), ", ", round(conf.high, 3), "]")
  )

b_value_gt %>%
  gt() %>%
  tab_header(
    title = "Estimated Allometric Slope",
    subtitle = "Mixed-effects model with tank as random intercept"
  ) %>%
  cols_label(
    Term = "",
    Estimate = "Estimate",
    `95% CI` = "95% Confidence Interval"
  ) %>%
  tab_options(
    table.font.size = 14,
    heading.title.font.size = 16,
    heading.subtitle.font.size = 13
  )

	Estimate	95% Confidence Interval
Estimated Allometric Slope
Mixed-effects model with tank as random intercept
Allometric slope (b-value)	0.99	[0.987, 0.994]

4. Surface Area Growth Analysis

Prepare SA-Adjusted Dataset

# Create dataset with surface area and time adjustments
bw_sa_df <- bw_sa_wide %>%
  drop_na(initial_weight, final_weight, wound) %>%
  filter(coral_id != 103) %>%  # Remove outlier
  mutate(
    delta_mass = final_weight - initial_weight,
    wound = factor(wound, levels = c("No Wound", "Small", "Large"))
  )

bw_sa_stats <- bw_sa_df %>%
  drop_na(initial_weight, final_weight, fish, wound, tank) %>%
  mutate(
    wound = factor(wound, levels = c("No Wound", "Small", "Large")),
    fish = factor(fish),
    tank = factor(tank),
    # Growth rate: mass change per surface area per day
    bw_sa_time = (final_weight - initial_weight) / (sa_cal * 21)
  )

cat("✓ SA-adjusted growth data prepared\n")

## ✓ SA-adjusted growth data prepared

cat("  - Valid observations:", nrow(bw_sa_stats), "corals\n")

##   - Valid observations: 66 corals

cat("  - Mean growth rate:", round(mean(bw_sa_stats$bw_sa_time, na.rm = TRUE), 5), "g/cm²/day\n")

##   - Mean growth rate: 0.00021 g/cm²/day

Visualize SA-Growth Relationship

sa_growth_plot <- ggplot(bw_sa_stats, aes(x = sa_cal, y = delta_mass, color = wound)) +
  geom_point(size = 2.2, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1.1) +
  scale_color_brewer(palette = "Dark2") +
  labs(
    title = "Coral Growth as a Function of Surface Area",
    subtitle = "Colored by Wound Treatment",
    x = "Surface Area (cm²)",
    y = "Mass Change (g)",
    color = "Wound Treatment"
  ) +
  theme_pubr(base_size = 14)

print(sa_growth_plot)

# Save figure
ggsave(
  filename = here("figures", "growth", "sa_growth_relationship.pdf"),
  plot = sa_growth_plot,
  width = 8, height = 6, dpi = 600, device = cairo_pdf
)

ggsave(
  filename = here("figures", "growth", "sa_growth_relationship.png"),
  plot = sa_growth_plot,
  width = 10, height = 7.5, dpi = 300, bg = "white"
)

cat("✓ SA-growth plot saved\n")

## ✓ SA-growth plot saved

5. Statistical Models: Allometrically-Scaled Growth

Prepare Scaled Growth Variable

# Extract b-value from allometric model
b_est <- b_value$estimate[1]

# Calculate allometrically-scaled growth
growth_stats <- bw_merged %>%
  drop_na(initial_weight, final_weight, fish, wound, tank) %>%
  mutate(
    wound = factor(wound, levels = c("No Wound", "Small", "Large")),
    fish = factor(fish),
    tank = factor(tank),
    # Size-corrected growth using estimated b-value
    growth_scaled = log(final_weight / initial_weight^b_est)
  )

cat("✓ Scaled growth variable calculated\n")

## ✓ Scaled growth variable calculated

cat("  - Using b-value:", round(b_est, 3), "\n")

##   - Using b-value: 0.99

cat("  - Mean scaled growth:", round(mean(growth_stats$growth_scaled, na.rm = TRUE), 4), "\n")

##   - Mean scaled growth: 0.0446

Fit Candidate Models

# Full model with interaction
mod_int <- lmer(growth_scaled ~ fish * wound + (1 | tank), data = growth_stats)

# Additive model (no interaction)
mod_add <- lmer(growth_scaled ~ fish + wound + (1 | tank), data = growth_stats)

# Fish-only model
mod_fish <- lmer(growth_scaled ~ fish + (1 | tank), data = growth_stats)

# Wound-only model
mod_wound <- lmer(growth_scaled ~ wound + (1 | tank), data = growth_stats)

# Null model
mod_null <- lmer(growth_scaled ~ 1 + (1 | tank), data = growth_stats)

cat("✓ All models fitted successfully\n")

## ✓ All models fitted successfully

Likelihood Ratio Tests

# Test for interaction
lrt_interaction <- anova(mod_add, mod_int, test = "Chisq")

# Test for fish effect
lrt_fish <- anova(mod_wound, mod_add, test = "Chisq")

# Test for wound effect
lrt_wound <- anova(mod_fish, mod_add, test = "Chisq")

# Compile results
lrt_table <- tibble(
  Test = c("Interaction (Fish × Wound)", "Fish Effect", "Wound Effect"),
  ChiSq = c(lrt_interaction$Chisq[2], lrt_fish$Chisq[2], lrt_wound$Chisq[2]),
  Df = c(lrt_interaction$Df[2], lrt_fish$Df[2], lrt_wound$Df[2]),
  p_value = c(lrt_interaction$`Pr(>Chisq)`[2], lrt_fish$`Pr(>Chisq)`[2], lrt_wound$`Pr(>Chisq)`[2])
) %>%
  mutate(across(where(is.numeric), round, 3))

# Format table
lrt_table %>%
  gt() %>%
  tab_header(
    title = "Likelihood Ratio Tests: Allometrically-Scaled Growth"
  ) %>%
  cols_label(
    Test = "Model Comparison",
    ChiSq = "χ²",
    Df = "df",
    p_value = "P-value"
  ) %>%
  fmt_number(columns = c(ChiSq, p_value), decimals = 3) %>%
  tab_options(
    table.font.size = 14,
    heading.title.font.size = 16
  )

Model Comparison	χ²	df	P-value
Likelihood Ratio Tests: Allometrically-Scaled Growth
Interaction (Fish × Wound)	4.225	2	0.121
Fish Effect	0.011	1	0.918
Wound Effect	1.856	2	0.395

6. Growth Rate Analysis (Mass/SA/Time)

Fit Growth Rate Models

# Full model with interaction
mod_int_sa <- lmer(bw_sa_time ~ fish * wound + (1 | tank), data = bw_sa_stats)

# Additive model
mod_add_sa <- lmer(bw_sa_time ~ fish + wound + (1 | tank), data = bw_sa_stats)

# Fish-only
mod_fish_sa <- lmer(bw_sa_time ~ fish + (1 | tank), data = bw_sa_stats)

# Wound-only
mod_wound_sa <- lmer(bw_sa_time ~ wound + (1 | tank), data = bw_sa_stats)

# Null
mod_null_sa <- lmer(bw_sa_time ~ 1 + (1 | tank), data = bw_sa_stats)

cat("✓ Growth rate models fitted\n")

## ✓ Growth rate models fitted

Growth Rate Model Comparison

# Likelihood ratio tests
lrt_interaction_sa <- anova(mod_add_sa, mod_int_sa, test = "Chisq")
lrt_fish_sa <- anova(mod_wound_sa, mod_add_sa, test = "Chisq")
lrt_wound_sa <- anova(mod_fish_sa, mod_add_sa, test = "Chisq")

# Compile results
lrt_table_sa <- tibble(
  Test = c("Interaction (Fish × Wound)", "Fish Effect", "Wound Effect"),
  ChiSq = c(lrt_interaction_sa$Chisq[2], lrt_fish_sa$Chisq[2], lrt_wound_sa$Chisq[2]),
  Df = c(lrt_interaction_sa$Df[2], lrt_fish_sa$Df[2], lrt_wound_sa$Df[2]),
  p_value = c(lrt_interaction_sa$`Pr(>Chisq)`[2], lrt_fish_sa$`Pr(>Chisq)`[2], lrt_wound_sa$`Pr(>Chisq)`[2])
) %>%
  mutate(across(where(is.numeric), round, 3))

# Format table
lrt_table_sa %>%
  gt() %>%
  tab_header(
    title = "Likelihood Ratio Tests: Growth Rate (g/cm²/day)"
  ) %>%
  cols_label(
    Test = "Model Comparison",
    ChiSq = "χ²",
    Df = "df",
    p_value = "P-value"
  ) %>%
  fmt_number(columns = c(ChiSq, p_value), decimals = 3) %>%
  tab_options(
    table.font.size = 14,
    heading.title.font.size = 16
  )

Model Comparison	χ²	df	P-value
Likelihood Ratio Tests: Growth Rate (g/cm²/day)
Interaction (Fish × Wound)	3.051	2	0.218
Fish Effect	0.058	1	0.809
Wound Effect	0.228	2	0.892

7. Publication Figure: Growth Rate

# Prepare data with correct factor ordering
bw_sa_stats$wound <- factor(bw_sa_stats$wound, levels = c("No Wound", "Small", "Large"))
bw_sa_stats$fish <- factor(bw_sa_stats$fish, levels = c("No Fish", "Fish"))

# Create publication-quality figure
growth_rate_plot <- bw_sa_stats %>%
  ggplot(aes(x = wound, y = bw_sa_time, color = fish, shape = fish)) +
  
  # Individual data points
  geom_jitter(
    aes(color = fish),
    position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.6),
    size = 2.5,
    alpha = 0.25,
    shape = 16
  ) +
  
  # Mean points
  stat_summary(
    fun = mean,
    geom = "point",
    size = 5,
    position = position_dodge(width = 0.6)
  ) +
  
  # Error bars (mean ± SD)
  stat_summary(
    fun.data = function(x) {
      data.frame(y = mean(x), ymin = mean(x) - sd(x), ymax = mean(x) + sd(x))
    },
    geom = "errorbar",
    width = 0,
    size = 1,
    position = position_dodge(width = 0.6)
  ) +
  
  # Scales
  scale_color_manual(
    values = fish_palette,
    name = "Fish Presence",
    labels = c("No Fish", "Fish")
  ) +
  scale_shape_manual(
    values = c(16, 18),
    name = "Fish Presence",
    labels = c("No Fish", "Fish")
  ) +
  
  # Theme
  theme_minimal(base_size = 16) +
  theme(
    panel.grid = element_blank(),
    axis.title = element_text(size = 20, face = "bold"),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 16),
    legend.position = "top",
    legend.title = element_text(size = 18, face = "bold"),
    legend.text = element_text(size = 16),
    legend.key.size = unit(1.5, "lines"),
    legend.margin = margin(0, 0, 10, 0),
    legend.box.spacing = unit(0.5, "lines"),
    axis.line.x = element_line(color = "black", size = 0.5),
    axis.line.y = element_line(color = "black", size = 0.5)
  ) +
  labs(
    y = "Growth Rate (mg/cm²/day)",
    x = "Wound Size"
  )

print(growth_rate_plot)

# Save figure
ggsave(
  filename = here("figures", "growth", "summary_line_growth_rate_by_wound_fish.pdf"),
  plot = growth_rate_plot,
  width = 7, height = 5, dpi = 600, device = cairo_pdf
)

ggsave(
  filename = here("figures", "growth", "summary_line_growth_rate_by_wound_fish.png"),
  plot = growth_rate_plot,
  width = 8, height = 6, dpi = 300, bg = "white"
)

cat("✓ Growth rate figure saved\n")

## ✓ Growth rate figure saved

8. Model Diagnostics

DHARMa Diagnostics: Growth Rate Model

library(DHARMa)

cat("=== Diagnostics for Growth Rate Model (bw_sa_time) ===\n\n")

## === Diagnostics for Growth Rate Model (bw_sa_time) ===

sim_growth_full <- simulateResiduals(fittedModel = mod_int_sa, n = 1000)

# Save diagnostic plot
png(
  here("figures", "diagnostics", "growth_dharma_full.png"),
  width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim_growth_full, main = "DHARMa Diagnostics: Growth Rate Full Model")
dev.off()

## quartz_off_screen 
##                 2

# Statistical tests
cat("\nDispersion test:\n")

## 
## Dispersion test:

print(testDispersion(sim_growth_full))

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91737, p-value = 0.67
## alternative hypothesis: two.sided

cat("\nOutlier test:\n")

## 
## Outlier test:

print(testOutliers(sim_growth_full))

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim_growth_full
## outliers at both margin(s) = 0, observations = 66, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.00000000 0.05435885
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0

Predicted vs Observed: Growth Model

bw_sa_stats$pred_growth <- predict(mod_int_sa)

p_growth_fit <- ggplot(bw_sa_stats, aes(x = pred_growth, y = bw_sa_time)) +
  geom_point(aes(color = fish), size = 2.5, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "darkred", linewidth = 1) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  scale_color_manual(values = fish_palette) +
  labs(
    x = "Predicted Growth Rate (mg/cm²/day)",
    y = "Observed Growth Rate (mg/cm²/day)",
    title = "Growth Model Fit: Predicted vs Observed",
    color = "Fish Presence"
  ) +
  theme_bw(base_size = 14) +
  theme(legend.position = "top")

print(p_growth_fit)

ggsave(
  here("figures", "diagnostics", "growth_fit_scatterplot.png"),
  p_growth_fit,
  width = 10, height = 8, dpi = 300, bg = "white"
)

DHARMa Diagnostics: Allometric Model

cat("\n\n=== Diagnostics for Allometric Scaling Model ===\n\n")

## 
## 
## === Diagnostics for Allometric Scaling Model ===

sim_allometric <- simulateResiduals(fittedModel = model_lmer, n = 1000)

png(
  here("figures", "diagnostics", "allometric_dharma.png"),
  width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim_allometric, main = "DHARMa Diagnostics: Allometric Scaling Model")
dev.off()

## quartz_off_screen 
##                 2

cat("\nDispersion test:\n")

## 
## Dispersion test:

print(testDispersion(sim_allometric))

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98124, p-value = 0.972
## alternative hypothesis: two.sided

cat("\nOutlier test:\n")

## 
## Outlier test:

print(testOutliers(sim_allometric))

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim_allometric
## outliers at both margin(s) = 1, observations = 67, p-value = 0.1254
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0003778063 0.0803764506
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.01492537

Model Summary Statistics

cat("\n=== Growth Model Diagnostics Summary ===\n")

## 
## === Growth Model Diagnostics Summary ===

cat("Residual standard deviation:", sigma(mod_int_sa), "\n")

## Residual standard deviation: 6.533866e-05

cat("Number of observations:", nobs(mod_int_sa), "\n")

## Number of observations: 66

cat("\nResidual summary:\n")

## 
## Residual summary:

print(summary(residuals(mod_int_sa, type = "pearson")))

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.674e-04 -3.278e-05 -1.770e-06  0.000e+00  3.712e-05  1.963e-04

9. Session Information

sessionInfo()

## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gt_1.0.0            ggpubr_0.6.1        broom.mixed_0.2.9.6
##  [4] DHARMa_0.4.7        lme4_1.1-37         Matrix_1.7-3       
##  [7] janitor_2.2.1       here_1.0.1          lubridate_1.9.4    
## [10] forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4        
## [13] purrr_1.1.0         readr_2.1.5         tidyr_1.3.1        
## [16] tibble_3.3.0        ggplot2_3.5.2       tidyverse_2.0.0    
## 
## loaded via a namespace (and not attached):
##  [1] Rdpack_2.6.4       rlang_1.1.6        magrittr_2.0.3     snakecase_0.11.1  
##  [5] furrr_0.3.1        compiler_4.5.1     mgcv_1.9-3         systemfonts_1.2.3 
##  [9] vctrs_0.6.5        pkgconfig_2.0.3    crayon_1.5.3       fastmap_1.2.0     
## [13] backports_1.5.0    labeling_0.4.3     utf8_1.2.6         promises_1.3.3    
## [17] rmarkdown_2.29     tzdb_0.5.0         nloptr_2.2.1       ragg_1.4.0        
## [21] bit_4.6.0          xfun_0.53          cachem_1.1.0       jsonlite_2.0.0    
## [25] later_1.4.3        broom_1.0.9        parallel_4.5.1     R6_2.6.1          
## [29] gap.datasets_0.0.6 bslib_0.9.0        stringi_1.8.7      qgam_2.0.0        
## [33] RColorBrewer_1.1-3 parallelly_1.45.1  car_3.1-3          boot_1.3-31       
## [37] jquerylib_0.1.4    Rcpp_1.1.0         iterators_1.0.14   knitr_1.50        
## [41] httpuv_1.6.16      splines_4.5.1      timechange_0.3.0   tidyselect_1.2.1  
## [45] abind_1.4-8        yaml_2.3.10        doParallel_1.0.17  codetools_0.2-20  
## [49] listenv_0.9.1      lattice_0.22-7     plyr_1.8.9         shiny_1.11.1      
## [53] withr_3.0.2        evaluate_1.0.4     future_1.67.0      xml2_1.4.0        
## [57] pillar_1.11.0      gap_1.6            carData_3.0-5      foreach_1.5.2     
## [61] reformulas_0.4.1   generics_0.1.4     vroom_1.6.5        rprojroot_2.1.0   
## [65] hms_1.1.3          scales_1.4.0       minqa_1.2.8        globals_0.18.0    
## [69] xtable_1.8-4       glue_1.8.0         tools_4.5.1        ggsignif_0.6.4    
## [73] grid_4.5.1         rbibutils_2.3      nlme_3.1-168       Formula_1.2-5     
## [77] cli_3.6.5          textshaping_1.0.1  gtable_0.3.6       rstatix_0.7.2     
## [81] sass_0.4.10        digest_0.6.37      farver_2.1.2       htmltools_0.5.8.1 
## [85] lifecycle_1.0.4    mime_0.13          bit64_4.6.0-1      MASS_7.3-65

Analysis Complete ✓

All figures and diagnostics saved to figures/ directory.

Coral Skeletal Growth Analysis: Allometric Scaling and Fish Effects

Hayden Vega & Adrian Stier

2025-10-28